Density functional theory of vortex lattice melting in layered superconductors: 

a mean-field— substrate approach 
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We study the melting of the pancake vortex lattice in a layered superconductor in the limit of 
vanishing Josephson coupling. Our approach combines the methodology of a recently proposed 
mean-field substrate model for such systems with the classical density functional theory of freez- 
ing. We derive a free-energy functional in terms of a scalar order-parameter profile and use it to 
derive a simple formula describing the temperature dependence of the melting field. Our theoret- 
ical predictions are in good agreement with simulation data. The theoretical framework proposed 
is thermodynamically consistent and thus capable of describing the negative magnetization jump 
obtained in experiments. Such consistency is demonstrated by showing the equivalence of our ex- 
pression for the density discontinuity at the transition with the corresponding Clausius-Clapeyron 
relation. 



I. INTRODUCTION 

The melting of the vortex lattice is one of the 
most striking aspects of the phenomenology of high- 
temperature superconductors^^. It is now both 
theoreticallji^SiL& and experimentally?^ established 
that enhanced thermal fluctuations render the Abrikosov 
flux-line lattice unstable to a vortex liquid over a large 
part of the B-T phase diagram in these materials. 
Much interest has been devoted to the melting transi- 
tion in extremely anisotropic materialsZiiiiiSii 3 . such as 
Bi2Sr2CaCu208. The large value of the anisotropy 14 
in these materials motivates their modeling in terms of 
a one dimensional array of magnetically-coupled two- 
dimensional superconducting layers with no inter-layer 
Josephson coupling. Vortex lines are then described by 
linear arrangements (stacks) of pancake vortices^, whose 
highly anisotropic interactions are mediated only by the 
magnetic field. 

The first theories for the melting of the vortex lattice 5 
were based on the Lindemann criterion^. The Lindc- 
mann melting criterion yields qualitatively accurate es- 
timates for the location of the first-order melting transi- 
tion. However, since it is based on the crystal properties 
only, it describes an instability rather than the melting 
itself. More detailed approaches, in particular analyses 
encompassing both solid and liquid phases, are thus re- 
quired for a quantitative description of the transition and 
of the discontinuities in thermodynamic quantities upon 
melting. 

Apart from numerical simulations^, the extreme 
limit of zero Josephson coupling between the layers 
has been analysed through classical Density Functional 
Theoriea 11 ! 18 ! 19 (DFT) and using a "mean-field" sub- 
strate approach 1 ^. Even if only the electromagnetic inter- 
action is retained, the analysis of the melting transition is 
a challenging task due to the long-range character of the 
interaction. At large magnetic fields the strong in-plane 
vortex repulsion dominates over the out-of-plane inter- 
actions and the melting line approaches the 2D melting 



temperature T^ D of the individual planes. On the other 
hand, at low fields, ignoring the possibility of low-field 
reentrance^, only few pancake vortex stacks are present in 
the system. Thermal fluctuations trigger the evaporation 
of isolated vortex stacks^ 5 - at T BKT , in correspondence 
to the zero-field two-dimensional Berezinskii-Kosterlitz- 
Thouless (BKT) transition. For intermediate values of 
the external magnetic field the melting line interpolates 
between the BKT transition temperature T BKT < T c 
close to the material critical temperature and the two- 
dimensional lattice melting temperature <C T c . In 
this regime, the full three-dimensional character of the 
system must be accounted for. The weak long-range out- 
of-plane interaction motivates a mean-field treatment*^, 
in which the 2D pancake vortex system in each layer is 
considered separately, with the effects of the other layers 
included via a self-consistent substrate potential. Exten- 
sive Monte Carlo and molecular dynamics simulations 
have confirmed the validity of this approximation and 
the accuracy of its results^. 

If the results of previous DFT approaches to the nu- 
merical melting line are compared with the simulation 
data, consistent results are obtained at magnetic fields 
larger than about Q.hB\ (B\ = $o/A 2 , where $0 = hc/2e 
is the unit flux and A the planar penetration depth). 
However, at lower fields, the disagreement between the 
two theories is substantial, with the DFT melting line 
shifted to much higher temperatures in comparison with 
the simulation data. Another difficulty with previous 
DFT studies is associated with the prediction of the sign 
of the magnetization jump across the transition. As it is 
well known, the vortex lattice exhibits a negative jump 
in magnetization upon freezing, leading to a solid phase 
which is less dense than the liquid, as in the ice- water 
transition. Such anomalous behaviour was first obtained 
in direct measurements of the magnetization discontinu- 
ity across the transition 9 . These observations are con- 
sistent with the negative slope of the melting line in the 
T-H phase diagram taken together with the Clausius- 
Clapeyron relation. However, earlier DFT analyses of 
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the freezing transition in the pancake vortex system ob- 
tained thermodynamically inconsistent results^, report- 
ing a positive density change and a negative slope of the 
melting line. 

In this paper, we adapt the DFT analysis to incor- 
porate ideas from the substrate approach; our analysis 
is simpler and more accurate in comparison with simula- 
tion data. Moreover, it addresses and solves the problems 
with the thermodynamic consistency faced in previous 
studies. The approach of this paper was first described 
in a recent letter where the effects of a surface on the 
melting transition have been studied^. In the present 
work we concentrate on the implications of our approach 
for the bulk transition in the infinite system. 

The main methodological difference between the work 
described here and previous DFT approaches lies in the 
determination of the direct correlation function. While 
in previous work, this correlation function was derived 
ab initio from the microscopic vortex interaction using 
the hypernetted chain approach or more elaborate exten- 
sions, here it is obtained by combining results from Monte 
Carlo simulations of 2D logarithmically interacting par- 
ticles, i.e., the One Component Plasma (OCP), with the 
substrate potential approach. The correlations of the 2D 
OCP are used to describe the effects of the strong in- 
plane vortex interaction, while the substrate potential 
accounts for the out-of-plane contributions. Within this 
new approach we obtain a simple expression for the free- 
energy which can be extended to the inhomogeneous case. 

We also prove the thermodynamic consistency of our 
approach by showing how to obtain the negative den- 
sity jump across the transition, fully consistent with the 
Clausius-Clapeyron equation. We do this by including 
a constraint which enforces an integer number of par- 
ticles per unit cell2i. Contrary to earlier claims^, no 
higher order correlation functions (three point or more) 
are needed to obtain the anomalous sign of the density 
discontinuity. 



II. MODEL 

Strongly anisotropic layered superconductors, such as 
the Bi- based compounds, are conveniently described in 
terms of the Lorentz-Doniach model^. The basic topo- 
logical objects are two dimensional vortices (pancake 
vortices) with a core limited to a single superconduct- 
ing layer. The interaction between pancake vortices is 
strongly anisotropic due to the underlying layered struc- 
ture. In Fourier space the potential readsi 



V(K, k z ) = 



<5> 2 d 2 



K 2 + k 2 



(1) 



4tt K 2 [l + \ 2 (K 2 + k 2 )}' 

here d is the layer spacing and A the bulk penetration 
depth in the plane. When placed on the same layer, 
pancake vortices feel a strong repulsive interaction, log- 
arithmically dependent on their separation, 



where £o = (<I>o/47rA) 2 is the vortex line energy and £ the 
(in-plane) correlation length. The interaction between 
vortices residing on different layers is attractive. Fourier 
transforming back along the z coordinate we obtain 



V Z {K)=- 



K 2 X 2 K, 



,-K+\z\ 



(3) 



with K + = ^J\/\ 2 + K 2 . Finally, transforming to planar 
real coordinates as well one obtains the potential in real 
space 



V Z (R) = ~ead 



+00 dR J (KR)e- K +M 



KK 4 



(4) 



For large in-plane separations R ^> A a logarith- 
mic attractive interaction is obtained V z ^q(R) ~ 
— (d/A)e~ z / A Vb(i?), suppressed by a factor d/\ when 
compared with the in-plane one. This out-of-plane in- 
teraction decays exponentially in the z direction over a 
distance A, extending over a large number (X/d) of layers. 
The strong anisotropy in the vortex interaction allows us 
to separate the strong in-plane repulsion from the weak 
out-of-plane interaction. The overall effect of the lat- 
ter can then be accounted for via an effective substrate 
termi^. 



III. CLASSICAL DENSITY FUNCTIONAL 
THEORY 

We follow the classical Density Functional Theory of 
freezing of Ramakrishnan and Yussouf described in Refs. 
I24II25I by choosing the uniform liquid as the reference 
state and estimating the difference in free energy due to 
the appearance of finite density modulations. For sim- 
plicity, we consider here a generic three-dimensional sys- 
tem of interacting particles; the modifications needed to 
describe the strongly anisotropic pancake vortex system 
are presented in the next section. 

The spatial arrangement of particles is described 
through the density field 



(r) = ^2S(r-n), 



(5) 



where is the position of the z-th particle (the index fj, 
emphasizes that p M (r) describes the non-averaged micro- 
scopic density) . To analyze the finite temperature behav- 
ior of the system we consider the density /?(r), averaged 
over thermal fluctuations 



P(r) = (p„(r)), 



(6) 



V (R) = -2e dMR/0, 



(2) 



where the brackets (. . .) denote the thermal average. The 
liquid and solid phases are characterized by qualitatively 
different density fields p(r): in the liquid phase the par- 
ticles are delocalized across the system and the averaged 
density p(r) = p 3D is constant; on the other hand the 
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solid phase is characterized by a modulated density p(r), 
with peaks at the lattice points. The appearance of finite 
density modulations is a consequence of particle-particle 
correlations arising from microscopic interactions. 

The classical DFT is based on the assumption that the 
free energy can be written as a functional of the averaged 
density p(r). One starts from the ideal gas free energy 
describing a non-interacting liquid and includes the cor- 
relations via an effective quadratic term in the density 
modulations Sp(r) = p(r) — p 3D . Within this approxima- 
tion the grand canonical free energy difference relative to 
the uniform liquid reads 



*"[p(r)] 
T 



d 6 r 



^ ld 3 r'8p(r)c(\r-r'\)8p(r') 



(7) 



where the temperature T is measured in unit of energy 
with fes = 1. The first two terms generalize^ the stan- 
dard free energy of an ideal gas to the case of non- 
homogeneous systems. The double integral term incor- 
porates the effects of interactions up to second order in 
the density difference Sp. This is the term which is re- 
sponsible for the appearance of finite density modulations 
which are absent in a non-interacting system. Therefore, 
the key input in this theory is the function c(r), the so- 
called direct pair correlation function, which accounts for 
correlations in the reference liquid. 

The direct pair correlation function can be related to 
more transparent physical quantities such as the static 
structure facto rial— i. Here, we give a brief derivation of 
the main relations which will be needed in our follow- 
ing discussion and, in particular, in the derivation of the 
Clausius-Clapeyron equation. We start from the micro- 
scopic density-density correlator 

(^( ri )^(r 2 )) = (Mr!)-p 3D ]Mr 2 )-^). (8) 

The Fourier transform of the density-density correlator 
defines the structure factor— 

5(q) = ^7 / d^d^e-^-^iSp^Sp^)), 

Next, we calculate the structure factor from the free 
energy J7J). The second functional derivative of the free 
energy 5il with respect to p(r) evaluated at p(r) = p 3D is 
the (functional) inverse of the density-density correlator 
(see also Ref. H|) 



6[6 P ( ri )]6[6p(r 2 )] 



= —S(r 2 -r 1 )-c(\r 2 ~r 1 \). (9) 

The functional inverse can be easily calculated in Fourier 
space. In this way, we find a relation between the struc- 
ture factor and the direct correlation function 29 



S(q) 



1 



1 - c(q) ' 



(10) 



Therefore, apart from additive constants (or delta func- 
tions in real space), the direct correlation function c(q) 
is given by the inverse of the structure factor S(q). The 
q = mode of the structure factor is related through the 
fluctuation-dissipation theorem to the isothermal com- 
pressibility of the system, i.e., S(q = 0) = ((-/V 2 ) — 
(N) 2 )/(N) 2 = p 3D Tn T (the isothermal compressibility 
of the ideal gas is = l/(p 3D T)) and thus 



1 - c(q = 0) = (p 3D TktY 



(11) 



The study of the free energy functional J7J) requires 
knowledge of the direct correlation function c(r) , a quan- 
tity which is usually obtained from the liquid state the- 
ory. It is possible to write a diagrammatic expansion 
for c(r) in terms of the microscopic two-body potential 
V[r) (or more precisely of the Mayer function f(r) — 
exp(— V(r)/T) — 1). For weak potentials, high-order cor- 
relations can be neglected and c(r) is given by the first 
(unperturbed) termSI 



c(r) « f(r) « -V(r)/T. 



(12) 



In order to obtain a better estimate for c(r), higher order 
terms in the perturbation expansion must be included. 
This is usually done by selecting specific classes of dia- 
grams out of the complete perturbative series. The ap- 
proach that is most widely used is the hypernetted chain 
(HNC) closure, an approximation scheme which, how- 
ever, is known to underestimate liquid-state correlations. 
The strategy we pursue here is different: we exploit the 
specific properties of the pancake vortex lattice by con- 
sidering the system as a collection of two-dimensional 
systems of log-interacting particles subject to a periodic 
modulated substrate potential due to the other layers, 
as done in Ref. 13. The function c(r) then combines an 
in-plane correlator, arising from the strong in-plane loga- 
rithmic interaction, and the weak out-of-plane potential. 



IV. DFT-SUBSTRATE APPROACH 

In contrast to standard liquids, the pancake vortex 
system exhibits a strong uniaxial anisotropy. Hence, we 
consider separately the in-plane (R) and out-of-plane (z) 
dependencies: in particular p(r) — > Pz(R) becomes a se- 
quence (in z) of two-dimensional densities. Similarly, we 
define the density variations Sp z (R) = p z (R) — p and the 
direct pair correlation function c z (R) (note that p is a 2D 
density). The DFT free energy of Eq. J7J can be adapted 
to the anisotropic vortex liquid 



Sn[ Pz (R)} _ [dz 2 
T Id 



p z (R) In EzQQ. - Sp z (R) 



1 [dz' 



d 2 R'^(R)c z _^(|R-R'|)^(R') .(13) 



The only input needed in the DFT free energy is the 
direct correlation function c z (R) which we obtain by im- 
plementing the substrate model for the pair correlation 
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function, 

Cz (R) = dc™(R)8(z) - (14) 

where V z (R) is the out-of-plane interaction of Eq. (0J ■ 

Within the planes, vortices are strongly correlated 
due to the repulsive logarithmic interactions (Eq. J5J). 
Hence, we can approximate cq(R) with the direct corre- 
lation function c 2D (R) of the two dimensional logarithmi- 
cally interacting particles (also known as one-component 
plasma, OCP). The in-plane component also contains 
contributions from the out-of-plane interactions V z=j ta (R) , 
which are small however, since V z ^o(R) appears in the 
perturbative expansion 2 ^ of cg(R) at least to quadratic 
order, [V Z ^ (R)/T] 2 ~ (d/X) 2 . Hence, the overall contri- 
bution due to out-of-plane interactions of all planes adds 
up to ~ (X/d)[V z ^(R)/T] 2 ~ (d/X). 

We use results of Monte Carlo simulations of the two- 
dimensional OCP at various coupling constants T = 
2sad/T to extract c 2D (i?). We have performed simula- 
tions on a system of 256 particles, using an alternative 
to the traditional Ewald summation method proposed 
recently by TyagiSi. Thermodynamic data and correla- 
tions were averaged over ~ 3 x 10 3 independent measure- 
ments following equilibration. The program was bench- 
marked using available numerical results for correlations 
and thermodynamic functions. In the simulations, the 
structure factor S(K) was calculated to yield the direct 
correlation function via the relation c(K) = 1 — 1/ 'S(K), 
see Eq. ifTUll . 

In the determination of the out-of-plane direct correla- 
tion function c z ^q(R) we neglect the higher order terms 
in the potential expansion, approximating it with the 
leading unperturbed value — V z ^q(R)/T (cf. Eq. fT2*|0 . 
Higher orders in c z ^q(R) involve at least terms of order 
~ (d/X) 2 . In the following, we will see that the relevant 
quantity in our analysis is the total out-of-plane corre- 
lator, defined as J(dz/d)c z ^o(R)- Whereas the leading 
term is of order (X/d)(d/ X) ~ 1 and hence comparable 
to the in-plane component, the subleading term is of or- 
der (X/d)(d/X) 2 ~ (d/X). We neglect this contribution 
to be consistent with our approximation for the in-plane 
component of the correlator. 

At a mean-field level the thermodynamically stable 
state corresponds to the minimal free energy configura- 
tion of the functional 1)13(1 . Then, the density functions 
p z (R) must obey the saddle point equation 

ln^ = J^-J d 2 R'c z ^ z ,(\R~R'\)SpAR'). (15) 

A key quantity in our discussion is the molecular 
field 24 ' 31 ! 32 £ Z (R) defined through 

e,(R)=Mp,(R)/p). (16) 

At the minimum of the free energy, combining the saddle 
point equation (|15|) with l|16fl . the molecular field be- 



comes 

6(R) =[^f[ d 2 R'<w(|R-R|)<W(R'). (17) 

Hence, the molecular field represents the screening poten- 
tial, produced by the modulated density which is acting 
back on the density itself. However, while 1(16(1 defines 
the molecular field and hence applies always, the inter- 
pretation in terms of a screening potential following from 
Eq. 1(17(1 is valid only at the minimum. 

V. FREE ENERGY IN FOURIER SPACE 

In thermodynamic equilibrium all superconducting 
planes are equivalent and the averaged vortex density 
p z (R) becomes independent of the layer position z; we 
write p z (R) = p(R)- Next, instead of seeking the exact 
form p(R) solving the non-linear integral equations 115(1 . 
we restrict our analysis to a simple family of periodic 
functions which model the modulations of the density 
in the triangular crystalized phase. In the following, we 
concentrate on the simplest case, retaining only the first 
Fourier components of the density in a triangular lattice 

^ = 1 + r, + V Me lKl - R = 1 + V + pg Kl (R), (18) 

p £r 

where the vectors Ki are the first reciprocal lattice vec- 
tors of the frozen structure and depend on the area a of 
the unit cell, p — Sp(Ki)/p is the Fourier component of 
the density with wave length K\, and rj — <5p(K = 0)/p is 
the relative density change upon freezing. Within this ap- 
proach the density is characterized by the three variables 
p, ?7, and K\. Note that, however, in a scheme where 
the number of particles in a unit cell is constrained to be 
an integer, the size of the unit cell a in the crystalized 
structure and the value of the density jump are related 
and, thus, Ki and r\ are not independent variables; we 
return to this issue in the next section. The function 

g Kl ( R ) = ^e lKl R = 2cos(2i) + 4cos(i) cos(y) (19) 

Ki 

includes the sum over the six first reciprocal vectors in the 
triangular lattice; in the last equality we have defined the 
dimensionless variables x = xKi/2 and y ~ \JiyK\/2. 
We also write a similar Ansatz for the molecular field 

£(R) = C + £WR), (20) 

retaining only the zeroth and first Fourier components, 
C and £ respectively, consistent with l(T7|l and the rapid 
decay of the correlator c(K) in the liquid phase, cf. Fig. 
□ 

The Fourier components of p(R) and £ (R) are not in- 
dependent and can be related through ((16|) . With the 
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FIG. 1: Direct correlation function at T/eod = 0.1 (T = 20) 
for the two-dimensional OCP, c 2D (_R") from MC simulations 
and c 2 ° b (if) (cf. Eq. 1(25(1 ) for the full three-dimensional pan- 
cake vortex system at the melting field B m (Bm/Bx ~ 0.099 
for T/eod = 0.1). We define our wave-length unit as G = 
(87r 2 p/v / 3) 1/ ' 2 • The substrate potential V sta <zk(K) modifies 
the full correlator c 2 ° b (K) in two different ways as compared 
to the in-plane correlation function c 2D (i("): i) it enhances 
the correlations of the liquid, pushing the melting to tem- 
peratures larger than T 2 ° and ii) it shifts the peak of the 
correlation function to a value of K which is smaller than 
G. In the incompressible limit K\ — G and at melting the 
full three-dimensional correlator assumes the critical value 
cJub = c 2 ° b (G) = c c « 0.856, cf. JUJl. Considering a finite 
compressibility, the system is allowed to gain correlation en- 
ergy by crystalizing at K\ < G, hence, at a density smaller 
than that of the liquid. However, a large density change is 
prevented by the finite compressibility of the system and the 
crystalized structure is characterized by a first reciprocal lat- 
tice vector Ki < G (see text). 



help of the relations 

- / d 2 Rg(K) = 0, - f d 2 R{g(H)} 2 = 6, (21) 

a Ja a J a 

we project the zeroth and first Fourier components of 
p(R) = pexp(C + and obtain 



C = -$(£)+m(l + 77), 

1 + 7? 



6 



■*'(£), 



where we have defined the function 



$(£) = In 



- / d 2 Re^( R ) 



(22) 



(23) 



This function is unaffected by rescaling of the unit cell 
area a (equivalently, it does not depend on K\). 

Substituting the Ansatze JTSJl and J5D|) in the DFT free 
energy (|13fl . we obtain the two-dimensional free energy 
density Suj 2 ° h = (d/pV) 5Vl as a function of the order pa- 
rameters r\ and p and the length K\ of the first reciprocal 



lattice vectors, 



(l+ry)[ln(l+ry)- $(£)]- 7? 
cL D b (0)77 2 



Hl(Ki)v, (24) 



where £ has to be understood as a function of r\ and p 
through (|22|l . In l|24|) we have defined the correlator 



dz V Z (K) 



T 



, (25) 



where V Z {K) — pV z {K) is the dimensionless Fourier 
transform (with an additional p factor) of the out-of- 
plane pancake vortex interaction (Eq. @). The out-of- 
plane interactions contribute to the correlator through 
the total stack potential 



T 



dzV z {K) 



Airpegd 



T TK 2 (X 2 K 2 + 1) ' 



(26) 



which enhances the nearest-neighbor correlations, cf. 

Fig. in 

Within this approach only the K = and the K = K\ 
components of the correlator are present in the expres- 
sion for the free energy H24fl . The K = Ki component 
measures the correlations which are responsible for the 
solidification of the liquid; its dependence on tempera- 
ture and field yields the position of the melting line and 
its local dependence on K determines the discontinuities 
in the first order transition. On the other hand, the 
K = component c 2 ° b (0) of the correlator is related to 
the compressibility of the vortex system via 1)11(1 . Using 
the relation kt = l/cn(0) between the compressibility 
and the compression modulus and using the expression^ 
c u (0) = B 2 /4n, we find that 



1 



,(0) = 



4:ne d B $ Bd 



T B, 



4ttT 



(27) 



consistent with the result of the present substrate-based 
approach: here, the K — ► divergence in the correla- 
tor c 2D (K — > 0) ~ —ATrpeod/TK 2 of the incompressible 
Coulomb gas^i is cancelled by the corresponding diver- 
gence in the out-of-plane component of the correlator, 
see (j2HJ, and the remaining term reproduces ((27(1 . In 
the Coulomb gas problem, this cancellation has to be 
achieved by the introduction of a compensating back- 
ground. 

Our expression for the free energy 1(24(1 in terms of 
the density p has to be compared with the original 
formula in Ref. I24L where the free energy has been 
given in terms of the molecular field £. The two ap- 
proaches are related via a Legendre transformation which 
adds the energy of an external periodic potential u(r), 
5W(u) = min p [SQ(p) — J d 3 ru(r)p(r)]. The relation 
((17(1 then includes the external potential u, £ Z (R) = 
u z (R) + (1/d) / dz' J d 2 R'c 2 _^(|R - K'\)5p z >(R'). Set- 
ting u = 0, both approaches provide identical results. 
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However, the formulation given here is more appropriate 
when describing non-uniform configurations as they oc- 
cur in a solid-liquid interface or near a surface, cf. Ref. 

M 



VI. CONSTRAINED THEORY 

Within our approximation for p(R), the state of the 
system is characterized by three parameters: the den- 
sity modulation p at the first reciprocal lattice vector, 
the density change 77 across the transition, and the wave 
number K\. However, a theory based on the Ansatz i|18|) 
and the functional l|13fl . in which 77, p, and K\ are inde- 
pendent variables, is not fully consistent. The problem, 
as was pointed out in Ref. |2l]j is that in l)18|l the density 
jump rj and the wave number K\ appear as independent 
variables. This theory therefore may lead to the appear- 
ance of states with a non-integer occupancy per unit cell 
and hence to an incorrect description of the solidification 
of the liquid 2 ^. To solve this inconsistency only states 
with a fixed integer total number of particles per unit 
cell should be considered. One must then proceed with 
a constrained minimization of the free energy, introduc- 
ing a Lagrange multiplier x which enforces the so-called 
'perfect crystal' conditional 



d 2 Rp(R) = 1, 



(28) 



i.e., each unit cell contains exactly one vortex. With this 
additional term the free energy density becomes 



6u,™(r,, p, K uX ) SwZ(r), P, Ki) 



T 
X_ 
pa 



d 2 Kp(R) - 1 , (29) 



where 5ui 2 ° b (r], p, K\) is given by i|24|). Substituting the 
Fourier Ansatz in the Lagrange multiplier term, we ob- 
tain the expression for the constrained free energy 



T 



X 



T 

(i + v) 



#i\ 2 
~G 



(30) 



where G = {^p/^/?,) 1 / 2 is the length of the first recip- 
rocal lattice vector associated with a solid with the same 
density as the liquid. 

To obtain the constrained minimum of the free energy, 
one needs to solve the saddle point equations for the vari- 
ables Xi Mi Vi an d K\. Taking the derivative with respect 
to the Lagrange multiplier we recover the 'perfect crys- 
tal' constraint (12811 in the form 



(31) 



i(*0 = (§) -'■ 



which yields a relation between 77 and K\. The case 
K\ = G describes an incompressible system, since the 
solid and the liquid have the same density. However, an 
ordinary first order phase transition is characterized by 
a finite jump of the density, and hence by a non-zero 77. 
Consistent with our constrained theory, a finite rj corre- 
sponds to the crystallization into a solid with a vortex 
density n v 7^ p and, hence, with a first reciprocal lat- 
tice vector K\ = (87r 2 n v /v / 3) 1 ^ 2 which is different from 
G, Ki 7^ G. When K\ > G the solid is denser than 
the liquid, as in conventional materials, whereas K\ < G 
leads to an anomalous density jump with a solid which is 
less dense than the liquid, as is the case in the water-ice 
transition. 

Next, the extremum condition of the free energy H30(l 
with respect to p gives the relation 



XKi)p, 



(32) 



which is essentially Eq. (|17|l in Fourier space. The equa- 
tion for rj 



x = -m + a-c™(o)H 



(33) 



is modified as compared with the standard unconstrained 
theory by the presence of the Lagrange multiplier, which 
can be found from the minimization of l|30|l as a function 
of K x 



X : 



3GV dc™ b {K) 



2Ki dK 



(34) 



Combining l|32|l - H34l) and the relation between p and f 
as given by (|22() . we can eliminate £ and x- The saddle 
point equations for rj and p can be written as 



M = 



[1 + V] 



f rf 2 R<?(R)exp[/ic 3 2 °(A' 1 )«7(R) 

J a 



1 



i-cL D b (0) 



d 2 Rexp A«£ub(*i)s(R-) 



[$(C b (^) M )- 



(35) 



3GV dc™ b {K) 



2K-, 



dK 



Ki 



(36) 



where 77 and Ki are related by (|31J) . The homogeneous 
(p = 0) and uncompressed (77 = 0) liquid always solves 
these equations. However, at low temperatures, other 
non-uniform solutions (p 7^ 0) may appear. If their cor- 
responding free energy is smaller than that of the liquid, 
5oj^° b < 0, the system freezes into the periodic (crystal) 
structure (note that at the minimum the perfect crystal 
constrained is fulfilled and 5Cj 2V , = Su) 2U , ). 

su b sub/ 

The equations l|35|) and H36fl can also be obtained di- 
rectly from the free energy l|24|) . 



8u™(i 1 ,p)^5u J ™ b {r!,p,K 1 {rj)), 



(37) 



where K\ is written as a function of 77 via l|31|) . The 
minimization of 5uj 2 ^ b (rj, p) /T with respect to 77 and p 
then yields the saddle point equations (|35|l and l|36|l . 
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0.04 




FIG. 2: Profiles of the two dimensional free energy difference 
5ui 2U (p) of Eq. 1381 as a function of the order parameter p, 
in correspondence to the values (from above to below) c 2D = 
c 2D (G) = 0.80, 0.83, 0.845, 0.856 (= c c critical, thicker line), 
0.87. At melting the order parameter jumps from the solid 
minimum at p s w 0.51 to p\ = 0, overcoming the barrier 
Suj™ w 0.0065 T. 



VII. INCOMPRESSIBLE LIMIT AND MELTING 

LINE 

The pancake vortex system is essentially incompress- 
ible (|c 2 ° b (0)| ^> 1) in a wide portion of the phase diagram 
and the corresponding value of the density jump is small, 
i.e., j) < 1. Hence, the determination of the melting line 
can be done within the incompressible limit, with rj = 
and K\ — G. This approximation is sufficiently accu- 
rate also for small magnetic fields (see later for estimates 
of the value of rj). In the incompressible limit, the free 
energy becomes a function of p alone, 

T T 

= 6£(/*)p - 3c 2 > 2 - §(£(/*)), (38) 

where c 2 ° b = c 2 ° b (G) and £(/x) is implicitly defined by 

m = $'(0/6- (39) 

The effective three-dimensional correlator c 2 ° h entering 
(|38[) is given by the sum of two contributions: the OCP 
correlation function c 2D = c 2D (G) and the stack potential 
V etaok (G)/T of Eq. JUJ, both evaluated at K = G. While 
the first depends on temperature only, the latter depends 
also on the vortex density and thus on the magnetic field, 



c™ b {T,B)=c 2D (T) + 
= c 2D (T) + 



'iirpsod 1 

TG 2 1 + A 2 G 2 
Vieod 1 



(40) 



2ttT [1+(8tt 2 /V3)B/B x ] 



where we use p/G 2 
B/B x . 



V3/{8ir 2 ) and X 2 G 2 = (8vr 2 /V3) 



At large values of B, the inter-plane interaction is 
negligible and the full correlator reduces to the 2D- 
OCP component c 2D . The temperature enters via the 
T-dependence of the direct correlation function c 2D (T), 
changing the coefficient of the quadratic term (as in the 
<j) A Landau theory). We obtain c 2D (T) directly from MC 
simulations; the results are shown in the inset of Fig. 
13 Increasing the temperature, the liquid correlations 
weaken, S{G) decreases, and so does c 2D , cf. Eq. (JTJJ. 

As a function of p, the free energy exhibits the shape of 
a Landau theory describing a first-order phase transition. 
In Fig. |21 we plot the free energy as a function of p for 
different values of c 2D . At large temperature, the correla- 
tor c 2D is small and 5u> 2V {p) exhibits only one minimum 
at p = with a value Suj 2D (0)/T = 0, in correspon- 
dence with the (homogeneous) liquid phase. Decreasing 
the temperature (which corresponds to increasing c 2D ), a 
second local minimum p s (metastable solid) with energy 



Suj 2D (p s ) , 2 



T 



0-2D A ,1^ ( -2D 

3c ^ s - $(c p s 



(41) 



appears in addition to the liquid minimum at p\ = 0. 
Freezing occurs when the liquid and solid minima assume 
the same value of the free energy, i.e., when Su> 2D (p s ) = 0. 
Within our single order parameter theory, this condition 
is equivalent to a simple equation for the correlator— 



0.856. 



(42) 



Going to even lower temperatures, c 2D further increases, 
the solid minimum decreases in value, 5ui 2D (p s )/T < 0, 
and the crystal becomes the only thermodynamically sta- 
ble phase. Monte Carlo simulations 34 show that the 
2D-OCP freezes at « Eod/70 where, however, the 
correlator assumes the value c 2D w 0.77 < c c (T 2 ° = 
2sod/T™ = 140). This disagreement is due to the ap- 
proximations we have adopted in our analysis. In partic- 
ular, at low temperatures, the higher-order peaks become 
important and more terms in the Fourier expansion have 
to be retained 26 . 

At lower magnetic fields, the inter-plane correlation 
becomes important and the 2D correlations c 2D are aug- 
mented by the stack potential V etack (G). The critical con- 
dition c 2 ° h (T,B) = c c can be solved together with 140(1 
and yields a simple expression for the melting line _B m (T) 



B m (T) 
B\ 



8tt 2 



V3r 



4tt(c c 



•CO) 



(43) 



This melting line is plotted in Fig. |3| (lower) together 
with the numerical results of the MC/MD simulations^. 
Furthermore, we find improved agreement in comparison 
with previous DFT studies where the direct correlation 
function was derived ab initio through approximative clo- 
sure schemes such as the hypernetted chain or the more 
elaborate Rogers- Young approachiiiS. In particular our 
novel approach approximates well the numerical results 
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FIG. 3: Comparison of the melting line obtained from the 
substrate-DFT analysis (solid line) with the result of full nu- 
merical simulations— (dashed line). Inset: values of the first 
peak at 7f ma x ~ G of the in-plane two-point direct correla- 
tor c 2D (K) as a function of T/eod = 2/T, from Monte Carlo 
simulations of the 2D one component plasma (OCP). 



for small values of the magnetic fields well (B < 0.5 B\), 
in a regime where results from previous work exhibit a 
substantial disagreement. 



VIII. DENSITY JUMP AND 
CLAUSIUS-CLAPEYRON RELATION 

To quantify the density jump across the transition, we 
need to consider the saddle point equation (|36[) for rj 



V 



*(£L(*i)/*) 

3G/i 2 dc™ b (K) 



2(1 + 77/2) dK 



Ki 



(44) 



where we have used K\ = GyT+ly » G(l+?7/2) to linear 
order in rj from i|31[l . Previous analyses^ of the freezing 
of the pancake vortex system were based on the uncon- 
strained free energy l|24l) . In the unconstrained theory 
the first reciprocal lattice vector is fixed at K\ = G and 
the equation for r\ contains only the first term in (|44p. 
since the second term is due to the Lagrange multiplier, 
cf. J2U). Hence, within the unconstrained theory, one 
obtains ('nc' stands for not-constrained) 



Vn 



iK&s) _ 4^r$(?> s ) 



i-C D b (o) 



<f>nBd 



(45) 



where we used O for c*°(0). Given that $(C>s) > 0, 
the sign of rj is always positive. This result contradicts 
the Clausius-Clapeyron relation and experimental evi- 
dence that vortices, like water, freeze into a solid which 
is less dense than the liquid. In a magnetic system the 
Clausius-Clapeyron relation reads^ 



AB 



-An As 



( dH m (T) yi 
\ dT ) ' 



(46) 



where AB = B\ — B s = —r\B\ = — ij&oP is the jump in 
magnetic induction and As = AS/V = (Si — S s )/V the 
(positive) jump in entropy density on heating. Ignoring 
the small difference between H and B, equation (|4(j[) can 
be rewritten as ('CC stands for Clausius-Clapeyron) 



An As (dB m (T) 
<S> p 



V dT J 



(47) 



Combining (|47|l with the negative slope of the melting 
line S m (T) of Eq. (|4*3l . we obtain that the density jump 
on heating is positive, which corresponds to a negative 
77, thus rjee < 0. Therefore, a theory with a positive rj 
and a melting line with a negative slope is not thermo- 
dynamically consistent. 

The second term in (|44H resolves this inconsistency. In 
2D systems the first maximum of the direct correlation 
function shows up at K max G, and hence dKC^ b {K)\c 
(and the second term in l|44|l ) is zero and rj > 0. However, 
for the 3D pancake vortex system the substrate poten- 
tial shifts the first peak of c™ b (K) to a value of K which 
is smaller than G, K max < G (cf . Fig. ^) . The system 
gains correlation energy by crystalizing at a K\ < G, 
hence, at a density which is smaller than that of the liq- 
uid. However, a large density change is prevented by 
the finite compressibility of the system and the crystal- 
ized structure is characterized by a first reciprocal lattice 
vector K\ below but still close to G, K\ < G. Hence, 
in (|44l) . the derivative dKC 2 J^ b (K)\K 1 at K\ is negative, 
d K c™ b {K)\ Kl « d K c™ b (K)\ G < 0, see Fig. □ The sec- 
ond term of l|44|) is therefore negative and a negative 
solution for 77 becomes possible. 

Next, we confirm the thermodynamic consistency of 
the constrained theory by showing that 14411 and the 
Clausius-Clapeyron equation (|47|l are equivalent. To 
compare (|47|l with (|44|l . we need to calculate the jump 
As in entropy density in l|47l) . The latter is given by the 
temperature derivative of the free energy difference 



P dSuj 2D b 
As = - - 3Ub 
d dT 



pT d Sco 2 ° b 



d dT T 



d dT 
d dT ' 



Tp dc 2 ° b (0) 2 
2d dT V 



(48) 



where we have used that Scu 2 ° b = along the melting line. 
The second term is of order rj 2 <C 1 and can be neglected 
when compared to the first one. In order to calculate 
the entropy jump and ?7cc, we need to evaluate the par- 
tial derivative dc 2 ° b (Ki)/dT at melting. The standard 
wa3 Tl8 i 31 is to estimate dc 2 ° b (Ki)/dT from the temper- 
ature dependence of the solid structure factor. Here we 
proceed differently. Comparing l|44|) with (|47l) (combined 
together with (|48|l ). we see that in 144|) the partial deriva- 
tive of c 2 ° b with respect to K\ appears whereas ijTT)) con- 
tains the partial derivative with respect to T, once we use 
(|48|l for the entropy jump. To compare the two different 
expressions for rj we need to find a way to connect these 
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two partial derivatives. This relation can be found from 
the critical condition which determines the melting line 
as we show in the following. 

The system freezes when the free energy at the solid 
minimum vanishes. Substituting the values of fj, s and of 
?7 S at the minimum into the expression for the free energy 
(|24[1 (or, equivalently, into (|30JO, the freezing condition 
reads 



T 



= 3C b (^)M s 2 ^ (1 + v s Wc™ b (KM = 0, 



where K\ is related to ?7 S through (|31|1 and we have ne- 
glected in H24f) the term (1 — c^ b (0))rj 2 , which is quadratic 
in t] s . At the minimum, the molecular field £ s and the 
order parameter /j s are related through £ s = c^ b (Ki)[x a , 
so the freezing condition can be rewritten as 



(3/cL D b (^i))4 2 -(i + ^)<i>(6) = 0. 



(49) 



For incompressible systems the same equation remains 
valid when one sets % = 



(3/cL D J4 2 - *(&) - o. 



(50) 



From the discussion in the last section we know that this 
equation is equivalent to the simple condition c*° b = c c , 
cf. 1(42(1 . Comparing (|50|l with (|49|l . it is easy to realize 
that the freezing equation in the compressible theory is 
equivalent to the one in the incompressible limit if one 
replaces cf° b (i"Ti)[l + rj(Ki)] by cj° b . Therefore, we can 
write a critical condition similar to 1(421) that is valid for 
a compressible system, 



(51) 



where we display explicitly both the T- and the in- 
dependences of the correlator. In 1)51(1 . the correlator 
depends indirectly on the magnetic field B in the solid 
phase through K\ via 



/8tt 2 B 
V\/3 $o 



1/2 



G(l+7?/2). 



(52) 



Consistently, since the magnetic field jumps across the 
transition (77 ^ 0), the magnetic field in the liquid B\ = 
<5>o/5 = B/(l + 77) is different from the magnetic field B 
in the solid phase. Along the melting line B m (T), K\ 
can be written as a function of T only by using l(52() . i.e. 
K x = K 1 (B nl (T)). Hence, at melting, the LHS of JSU 
can be written as a function of the temperature alone; 
taking the total derivative d/dT of l|51|l with respect to 
the temperature, we obtain 



dc™ b (T, Ki) _ c™(T, Kx) dr,{K x ) dK x dB n 



dT 



1 + 77(^1) OKi dB dT 
dcZiT^KjdK, dB m 



dK x dB dT 

We need to compute the derivative 

dKx _Ki_ G 

~dB ~ 2B ~ 2Si(l + r//2)' 



(53) 



(54) 
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FIG. 4: Values of the density jump 77 across the transition 
as a function of temperature for the 3D vortex system. The 
density jump 77 is negative and small, of order 10 -4 at low T 
(large magnetic fields) rising to 10 -2 at larger T (low magnetic 
fields). Inset: log-plot of the absolute value of r\. 



where we have used B = B\(l + 77) and the linearized 
relation between K\ and G in <|52|l . Inserting Eq. I|53|) 
(with the help of JH3J and lj3"T)l j into we obtain the 
entropy jump across the transition 



As = 



3/j 2 pT 
dB\ 



XT,Kx) 



1 



G 



2(1 + 77/2) OK 



dB m 
dT 



(55) 



Inserting this expression into l|47ll , we obtain the density 
jump 77CC described by the Clausius-Clapeyron relation 



1 



Vcc 



+ 



cL D b (0) 

3Gfi 2 



3c^ b (T,K 1 )/i 2 
I + 77 
dc^{T,K) 



2(1 + 77/2) 



dK 



(56) 



where in the last line we have used dBf/AirpT = 
d$oBi/47rT = 1 - cj° b (0) from QZJ. The first term in 
the square brackets in (|56|l can be rewritten with the 
help of the freezing condition. As a final result we obtain 



1 



Vcc 



,,,,, HcZ(Ki)v) 
3G/J 2 dc 2 ° b {T,K) 



2(1 + 77/2) dK 



V, 



(57) 



which is exactly Eq. (|44|1 . Thus, we conclude that the 
saddle point equation l|44|l is fully consistent with the 
Clausius-Clapeyron relation. 

To obtain an estimate for the density jump across the 
transition we have performed a numerical minimization 
of the constrained free energy (|30|l . The system freezes 
when the solid minimum exhibits the same free energy 
as the liquid phase, i.e., when 6uj(r] s , fj, s ) — 0. In Fig. 
0] we present the results of our numerical analysis for 
various values of T. For each T we show the density 
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jump at the transition; as expected, we find a negative 
value of 77. However, the modulus of the density jump rj 
is always small, between \<q\ « 10~ 4 at low temperatures 
(large B) and \r)\ ~ 10 -2 for large temperatures (low B). 
The effect of such a small r\ on the determination of the 
melting line and on the value of \x s in the solid phase 
is negligible. Finally, from l|47|l we can obtain the value 
of the entropy jump across the transition. In agreement 
with previous studies^, we obtain a density jump per 
pancake vortex AS/N ss 0.4 fc B at large fields decreasing 
to AS/N w 0.1 /cb at low magnetic fields corresponding 
to temperatures of order T = 0.2 eqCL. 

IX. CONCLUSIONS 

We have presented an analysis of the freezing transi- 
tion of the magnetically coupled pancake vortex system 
within a classical density functional theory. Despite the 
simplicity of our approach, our results represent a con- 
siderable improvement when compared to the predictions 
of earlier work, particularly in the determination of the 



melting line. Moreover, we have addressed the prob- 
lems with the thermodynamic inconsistency which af- 
fected earlier work. We showed how to obtain a negative 
density jump in accordance with the Clausius-Clapeyron 
equation and the retrograded melting line. Note that our 
derivation of the Clausius-Clapeyron equation is not lim- 
ited to the present vortex system but can be generalized 
to non-magnetic systems. 

The techniques described in this paper are easily ex- 
tended to the study of inhomogeneous situations, as was 
already done in the analysis of surface effects on the melt- 
ing transition 2 ^. The study of artificial surface pinning 
potentials offers itself as another application of our DFT 
approach. Such an analysis would shed light on the re- 
sults of recent experiments described in Ref. |35| where 
the response of vortices in BiSCCO to a weak perturba- 
tion induced by pinning structures created on the sample 
surface has been investigated. 

We acknowledge support from the Swiss National 
Foundation through MaNEP (ADC) and from the DST 
(India) through a Swarnajayanti Fellowship (GIM). 
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